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Abstract 

We construct an extension of the standard Kogut-Susskind lattice model for classical 3+1 
dimensional Yang-Mills theory, in which "classical particle" degrees of freedom are added. 
We argue that this will correctly reproduce the "hard thermal loop" effects of hard degrees 
of freedom, while giving a local implementation which is numerically tractable. We prove 
that the extended system is Hamiltonian and has the same thermodynamics as dimensionally 
reduced hot Yang-Mills theory put on a lattice. We present a numerical update algorithm 
and study the abelian theory to verify that the classical gauge theory self-energy is correctly 
modified. Then we use the extended system to study the diffusion constant for Chern-Simons 
number. We verify the Arnold-Son- Yaffe picture that the diffusion constant is inversely 
proportional to hard thermal loop strength. Our numbers correspond to a diffusion constant 
of T = 29 ± 6c^T 4 for ml = llg 2 T 2 /Q. 

PACS numbers: ll.10.Wx, 11.15.Ha, ll.15.Kc 

1 Introduction 

One of the most profound and most poorly explained observations of modern cosmology is 
that the universe contains macroscopic amounts of matter, but not of antimatter. Since 
baryon number is conserved in all observed terrestrial and astrophysical phenomena, the 
existence of such a baryon number asymmetry (baryons over antibaryons) seems peculiar, 
particularly because it clearly violates C and CP symmetry. Further, the abundance of 
baryons, as compared to the abundance of photons or the entropy density of the universe, 
is a remarkably small number, on order 10 -10 ; and since the entropy of the universe has 
changed very little since its very hot early epochs, this small nonzero number is an initial 
condition to the universe viewed at least back to the epoch of primordial nucleosynthesis. 

Sakharov [|IJ made the first attempt to understand what is involved in explaining this 
observation. If the universe does not begin with such an asymmetry, then baryon number 
must be violated to generate it; so must C and CP. He also noted that if baryon number 
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is not conserved, then in equilibrium it will go to zero, as a consequence of CPT symmetry. 
Hence, if the baryon asymmetry of the universe was generated dynamically, the universe must 
in its early history have gone through some departure from thermal equilibrium, in which 
baryon number violation was active and C and CP violating physics were relevant, followed 
by a steep suppression of the rate of baryon number violation before thermal equilibrium 
resumed (to prevent the baryon number from being erased again). And if baryon number 
was violated at a rate faster than the Hubble expansion at any point in the universe's history, 
then the baryon number must have been generated dynamically. 

While it is known that grand unified theories (GUTs) generically violate baryon number, 
it turns out that baryon number is violated already in the standard model, as first shown by 
t'Hooft 0. The key observation is that baryon number is coupled through the axial anomaly 
to the Chern-Simons number (Nqs) of the SU(2) weak hypercharge field, 

^ = 7V F ^, (JV F = 3), (1) 

and that iV cs change is thermally activated in the hot electroweak plasma @, f|, ||. This 
baryon number violation occurs at much lower temperature ranges, and hence later in the 
evolution of the early universe, than GUT mechanisms; it could therefore erase any baryon 
number asymmetry generated when GUT mechanisms are relevant^ and generate the ob- 
served abundance during the electroweak phase transition, when the rate abruptly shuts off. 
This latter possibility has stimulated the field of electroweak baryogenesis. 

One thing we need to know to understand electroweak baryogenesis is how quickly baryon 
number is violated in the standard model, above the electroweak phase transition. The 
violation involves nonperturbative physics, and there are no known analytic methods which 
are reliable above the phase transition temperature (in contrast to the situation below it 
0, Q). Already at the thermodynamic level the fields responsible for the violation of baryon 
number suffer from the "infrared problem" of thermal field theories with light interacting 
bosons. 

The "infrared problem" of the thermodynamics of thermal field theories with light in- 
teracting bosons can be "solved" by the dimensional reduction procedure || [7|, ||, which 
reduces the problem to a three dimensional path integral which can be put on a lattice 
H [H| |ll|] and studied numerically (TJ, [TJ, [L4], |TJ, [TjJ. This "dimensional reduction" turns 



out to be identical to considering the partition function of the 3+1 dimensional, classical 
bosonic theory, with certain counterterms [IjJ. The most hopeful approach to the study of 



infrared sensitive, dynamical properties of the plasma is that they also behave essentially 



classically |[Lq| , and there is reason to believe that they do |L5], [19 . 

There is a problem with the classical approximation to the dynamics, however. The 
ultraviolet modes (p ~ irT) certainly do not behave as classical fields, and they interact with 
the infrared modes (which we will assume do behave essentially classically). The classical 
approximation will serve only if the interaction between soft and hard modes is somehow 
unimportant to the evolution of the soft modes. It is known that, at weak coupling, the 
hard modes only influence the thermodynamics of the soft modes by shifting the Higgs mass 

1 The standard model violates B + L but not B — L. Extended GUTs such as SO(10) can produce a B — L 
excess, which electroweak processes will not touch. But the possibility of B production in the standard 
model certainly motivates the study of these scenarios. 
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and the Debye screening mass for the Aq field. f] However, the unequal time generalization of 
Debye screening, the hard thermal loops, constitute very nontrivial physics, physics which the 
lattice implementation of the classical theory gets wrong p0| , both in the size of the screening 
effects (which depend inversely on the lattice spacing) and in the functional details. 

Does this problem matter to the rate at which baryon number is violated? In particular, 
since there is a fluctuation dissipation relation which relates the rate at which a baryon 
number excess is dissipated in the plasma to the rate of Nqs diffusion per unit volume JH], 
we might ask whether the hard thermal loops matter to the Nqs diffusion constant. The 
numerical lattice results of Ambj0rn and Krasnitz implied that there is a well defined small 
lattice spacing limit to the classical, lattice diffusion constant |17|, but the definition of lattice 
Ncs used there gives bizarre results in the broken phase [^] due to serious lattice artifacts 



15]. Theoretical arguments suggest that hard thermal loops should slow down the evolution 



of infrared magnetic fields, for the same reason that infrared magnetic fields in an abelian 



plasma become pinned by the conductivity of the plasma [23, 24 ; hence one should expect 



that the Ncs diffusion constant vanishes linearly with lattice spacing in classical simulations. 



Recently, numerical techniques based on a topological definition of iV cs |25[ and a "cooled 
field" definition which removes the worst problems of the old definition |26 have indicated 
a lattice spacing dependence in the diffusion constant, though it does not appear to be as 
large as expected. Hence, at present the situation seems confused, and it certainly seems 
necessary to find a way to include more faithfully the effects of hard thermal loops in the 
study of iVcs diffusion, and in other studies of infrared bosonic field evolution, such as the 
study of the phase transition dynamics IJT5J and the generation of quark gluon plasma in 
heavy ion collisions [f27[ pj9fl . 

In this paper we develop, implement, and study a proposal by two of us J50] for including 
the hard thermal loop effects into classical, lattice simulations of Yang-Mills or Yang-Mills 
Higgs theory. The implementation produces a local, Hamiltonian system in which added 
"particle" degrees of freedom convey the hard thermal loop effects. We explicitly verify that 
the enlarged system conserves energy and symplectic measure, has the same thermodynamics 
as the quantum theory in the dimensional reduction approximation, and produces the right 
dynamical behavior for infrared fields. We examine in detail the retarded propagator of the 
abelian theory with particles to show that they correctly produce all of the physics of hard 
thermal loops, including Debye screening, plasma oscillations, and Landau damping. Then 
we apply the particle method to determine the Ncs diffusion constant and its dependence 
on the Debye mass. The results verify the arguments of Arnold, Son, and Yaffe (ASY). In 
particular, we verify that as the total strength of hard thermal loop effects varies by a factor 
of 3, with the lattice spacing and the physical volume held constant, the diffusion constant 
for iVcs also varies by a factor of 3. 

The outline of the paper is as follows. In Section [2] we review what the hard thermal loops 
are and how they can be understood in terms of Vlasov equations, ie in terms of the influence 
of a bath of particles on the classical infrared modes. In section we explicitly construct a 
spatial lattice, continuum time system which can be viewed as an N body simulation of the 
Vlasov equations. We show that this system is energy conserving and Hamiltonian and argue 



2 There are also perturbatively computable corrections to the couplings of the theory, which amount to a 
specification of the renormalization point, see Farakos et. al. ^ B 
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that, for small charge and large particle number, it reduces to the (lattice) Vlasov equations 
in the eikonal approximation, which are known to produce the right hard thermal loops. In 
Section |] we study its thermodynamics, which are shown to be the same as dimensionally 
reduced Yang-Mills theory put on the lattice, at a specific Debye mass set by the number 
and charge of the particles. Section |5| develops a stable, time symmetric update algorithm, 
and Section || numerically studies the retarded propagator of the abelian theory. Section 
|7] uses the evolution algorithm, together with the topological definition of N C s developed 
in [p5| , to study the N cs diffusion constant and its dependence on lattice spacing and on 
the size of hard thermal loop effects. Section |8] concludes. There are also four technical 
appendices. Appendix A analytically studies how the lattice nature of the electric fields 
changes the plasma frequency, in the limit of many particles of small charge. Appendix B 
studies the thermodynamics of the lattice system when ~ g 2 T/a. Appendix C reviews 
how to probe the retarded propagator via linear response to an external current. Appendix 
D discusses the relation between lattice and continuum time scales. 



2 Hard thermal loops and particles 

In this section we briefly review an approach to hard thermal loops due to Kelly et al |3l| and 



discuss how one can generate the influence of hard thermal loops by coupling the infrared 
fields to a heat bath of particles. 

In the diagrammatic approach to field theory at high temperature it was shown by 
Braaten and Pisarski |32| that a resummation procedure is needed in order to take into 
account consistently all contributions at leading order in the coupling constant. Such con- 
tributions are called "hard thermal loops" since they arise from one-loop diagrams with soft 
external legs and hard internal momenta. The hard thermal loop resummation was shown to 
produce gauge invariant results for physical quantities such as the gluon damping rate in a 
QCD plasma p3 |. An effective action for the hard thermal loops was derived by Taylor and 



Wong P|[] by imposing gauge invariance on the generating functional. Another approach 



developed by Blaizot and Iancu [|3f| is based on a truncation of the hierarchy of Schwinger- 
Dyson equations and the generating functional was obtained through a consistent expansion 
in the coupling constant. 

While remarkably insightful, the approaches mentioned above are quite technical, often 
involving lengthy calculations, and they tend to hide the classical nature of hard thermal 
loops. The hard thermal loops arise from loop diagrams and can be obtained from the 
Schwinger-Dyson equations of quantum field theory. But they are also UV finite, with 
loop integrals being exponentially suppressed in the ultraviolet. This is because they arise 
entirely from thermal fluctuations. One usually thinks of such fluctuations, at least in the 



ultraviolet, as being well described by classical particles, which motivated Kelly et al J31 
to find an alternative, classical formalism for hard thermal loops, by considering the linear 
response of an ensemble of thermal particles to a weakly coupled, slowly varying mean field in 
the framework of classical transport theory [j3f|. They start by considering particles carrying 



non-abelian SU(iV) charge q a , a = 1, N 2 — 1. The Wong equations |37| are used to describe 
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the proper time evolution of a particle with phase space coordinates (^,p' Jl ,q a ): 



m 



ni 



m- 



d^ 
dr 
dp* 1 

dr 
dq a 



P* 1 , 



gq a F^p u , 



(2) 
(3) 
(4) 



where g is the coupling constant and F£ v denotes the strength of the mean color field A%. 
Note that the color charge is itself subject to dynamical evolution, a feature which is absent 
in the abelian case. Both the dynamical evolution of a particle's spin and the spin coupling 
are neglected since spin interactions are down by an order of g in a weak mean field describing 
soft excitations (k ~ gT, ~ T, F^ u ~ gT 2 ) [[3^ . 

Consider the classical one particle distribution function f(x,p,q) evolving in time ac- 
cording to the Boltzmann transport equation, 



df(x,p,q) 

m : = U 

dr 



(5) 



where the collision integral describing hard collisions between particles is neglected since soft 
collisions mediated by the mean field dominate. Inserting the Wong equations (0)-(|j) into 
(^), one arrives at the following: 



V 



d 
dx^ 



gq a F? 



> w dp v 9UbcA » q dq- 



f(x,p,q) = 



(6) 



A self-consistent set of non-abelian Vlasov equations for the distribution function and the 
mean color field can be obtained by augmenting (Qj with the Yang-Mills equations, 



D,F^ = gj dpdqp»qf(x,p, q) = f(x) 



(7) 



The non-abelian Vlasov equations @ and (0) are now applied to study the soft excitations 
in a hot, color-neutral plasma with massless particles. In the spirit of linear response theory, 
one expands the distribution function in powers of g, 



f = r> + gf w + 0{g 2 ) , (8) 

where f^(po) = Cn BF (p ) is the equilibrium distribution in the absence of a net color field. 
At leading order in g, the equation (Q) reduces to 



V 



9 f A b c d 

gjabcA^q 7^ 



dx^ 



d 



f W (x,P,q)=p"q a F^—f^(p ) . 



Similarly, there is a net induced current density in momentum space, 

f a (x,p)=g 2 fdqp>>(ffW(x,p,q) . 



(9) 



(10) 
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From (H) and (|10D, a constraint satisfied by the color current can be derived. Using the 
standard field theory definition of effective action j^{x) = — /J^^ , one then arrives at the 
hard thermal loop effective action of the following form: 

rIO 1 

(11) 



r 1 2 

t H tl = 



(fxA a Q (x)A a Q (x) - / -^-W(A ■ v) 



(2nY 

where v = (l,p/p ), m D = gT J '^E± n ^+ 2N f j5] anc i integration / dfl is over all directions 
of the unit vector p/po. The first term describes Debye screening. W(A ■ v) in the second 



term is a functional. Its explicit form has been given by Taylor and Wong pi |, also by 



Efraty and Nair f39|. Note that the derivation of Thtl stays completely within the classical 



transport theory and makes no use of the quantum theory. This justifies the statement that 
hard thermal loop effects are classical (in the sense of the classical particle approximation). 



Also note that the form of the expression inside the brackets in Eq. (11) does not rely on the 
particles obeying any specific statistics; they can be Fermi, Bose, or even Boltzmann particles. 
Nor does it depend on what group representation they are in. These only affect how much 
a species contributes to the leading coefficient m^, which is given by a sum over charged 
species of their individual contributions. The ratio of a particle species' contribution to m 2 ^ 
and the mean density of those particles is 2g 2 C 2 (R)(E~ 1 ). (For classical, distinguishable 
particles in the adjoint representation, g 2 C 2 (R) is replaced by the mean squared value of the 
particle charge q in one Lie algebra direction.) 

It is worth noting that the magnetic field does not play a role in producing hard thermal 
loops at leading order. This can be easily seen from (^|), where only the electric field enters the 
term on the right. Physically, this is because the magnetic field just rotates the momentum 
distribution, and at leading order that distribution is the rotationally invariant thermal 
distribution. The magnetic field only influences existing departures from equilibrium, which 
is a subleading effect. In contrast, the electric field can polarize the plasma and generate a 
net current, which in turn interacts with the mean field and hence generates the desired hard 
thermal loops. Therefore, in solving the Wong equations, one can leave out the magnetic 
term in the Lorentz force in (|^) if the plasma is only slightly driven out of equilibrium. 
However, if one is interested in dynamical processes occurring in out-of-equilibrium plasma, 
then the force due to the magnetic field has to be included. 



The effective action ( JlT]) is conceptually simple and formally appealing. It provides a 
concise way of summarizing hard thermal loops and allows one to better understand the 
influence of the hard thermal modes on the soft excitations. Nevertheless, it does not prove 
a ready starting point for practical calculations. In particular, its nonlocality makes it hard 
to apply to study nonperturbative physics such as Chern-Simons number diffusion, where 
analytical methods are rare and one has to rely on numerical simulations. However, the 
fact that the hard thermal loop effective action ( |Tl| ) can be derived from classical transport 
theory, i.e. the Vlasov equations, implies that hard thermal loops can be generated by 
solving the Vlasov equations numerically. One could do so by an N body simulation by 
solving the coupled system of Wong equations (@)-(|) and Yang-Mills equations (|7|). This 

3 Here, N s is the number of fundamental scalars, and as in the preceding text, Np is the number of 
generations, each containing four chiral doublets; in QCD 2JVf would be replaced by Np the number of 
flavors, and N s would only appear in supersymmetric extensions. 
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has the advantage of being local in spacetime and allows a practical, real time study of both 
equilibrium and nonequilibrium phenomena. 

Based on these ideas, our goal will be to derive a system of dynamical equations which, 
in its low-frequency, long- wavelength limit, reduces to the same effective field theory as the 
full thermal gauge theory. We have no intention of modeling the gauge theory as truly as 
possible in all respects; rather, we will be satified to reproduce its infrared behavior correctly. 
For example, it will be irrelevant for us whether the system of classical particles carrying 
gauge charge has detailed resemblance to the hard thermal modes of the gauge field, if only 
its back- reaction produces the nonlocal HTL effective action (|TT|). The system of classical 
dynamical equations to be derived in the following sections will have three parameters: the 
classical thermal length scale (^ 2 T) _1 , the gauge charge (gQ) 2 of the particles, and the 
particle density (n). It will be sufficient for our purpose if we tune these parameters so that 
they produce the correct values for the parameters in the HTL action. 

3 The lattice system 

To reproduce the hard thermal loop effects, our idea is then to numerically implement "par- 
ticles" obeying Wong's equations and moving in the background of classical, lattice fields. 
Eventually we want a discrete time update algorithm, but an important step is to construct 
a continuous time, lattice system with the right properties. The general philosophy of adding 
particles with adjoint (Lie algebra) charge q of fixed magnitude, kinetic momentum p, and 
continuous position coordinate £ satisfying ^ = Pi/\p\, has been presented in [pD[ ; here we 
will specify the complete implementation. Our requirements for this system are 

• There must be conservation of energy, 

• The evolution should preserve the phase space (symplectic) measure, 

• The system should respect cubic, (lattice) translation, C, P, and T symmetry, 

• Gauss' law should be identically preserved, 

• The small lattice spacing a limit (or smooth field, large time limit if one thinks in 
lattice units) must recover the Yang-Mills field equations supplemented with Wong's 
equations, 

• The thermodynamics of the infrared classical fields must be given by the path integral of 
quantum Yang-Mills theory in the dimensional reduction approximation (or its lattice 
discretization). 

Note that the first two conditions ensure that the system is Hamiltonian, and hence that 
the thermodynamics are well defined. We would also like a thermalization algorithm for the 
system. 

We emphasize that the particles are a device to reproduce the hard thermal loops and 
should not be taken literally as reproducing all the behaviors of the hard modes. In particular 
it is not a problem that they are distinguishable, and that their number is conserved. We 
expect them to satisfy Boltzmann statistics, rather than Bose or Fermi-Dirac statistics, but 
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this is also not important to reproducing the functional form of the hard thermal loops, as 
discussed in the last section; we need only make sure that the number density and charge 
of the particles yields the desired value for m^. However, any hope that the method can be 
enhanced to account for physics beyond hard thermal loops is clearly remote. Hence we will 
not be able to say anything about effects which are subleading in g. 

Let us list the degrees of freedom of the proposed system (see Fig. |l] for illustration). 
We consider a 3-torus of spatial extent N 3 (in lattice units, which will be used throughout 
except when it is convenient to write the lattice spacing explicitly). "The lattice" will refer 
to the integer lattice on this space, ie all points with all three coordinates an integer. A link 
x, i will refer to the line between the lattice site x and the site x + i (henceforward x + i) and 
on each link there will be a parallel transporter Ui(x) G SU(2) and an electric field Ei(x) G 
LSU(2), the Lie algebra of SU(2). For each index value a G {l,..,N p } there is a particle 
with coordinate £ a defined on the torus, momentum p a G 3? 3 , and charge q a G LSU(2) 
satisfying q 2 a = Q 2 . By definition Ei(x) will be the left acting covariant time derivative of 
Ui(x), D Ui(x) = E^(x)iT a Ui(x), and the momentum p a will tell the direction the particle 
moves in, = p a ,i/\p\ a (note the are defined mod N). It remains to define update rules 
for E, p, and q. 

First, what is the meaning of q? It should be a charge which the classical fields U, E 
"see," but such a charge should reside at a lattice point. We take q a to "live" at the lattice 
point Xi closest to £ a , ie \xi — £ a ,i\ < 0.5. It will gauge transform as an adjoint object at that 
site; and the charge observed by the classical fields will be 

Kx) = V? a x( I ^ at f + . (12) 
a 1 £ not at x v ' 

This is the quantity which enters Gauss' law, 

p(x) = Ei{x) -Uj{x- i)E i {x - i)Ui{x -i) = D h - E{x) . (13) 

The current jf, an adjoint vector field, should be defined on the links and is specified by 
the requirement that it form a conserved current, 

p = U K X ~ ~ i)Ui{x - i) - ji(x) . (14) 

i 

The easiest way to see how this fixes the current is to think of the charge q a of particle a 
as living at site x until the particle moves to be closest to another site, say x + i; then the 
charge must abruptly slide along the link connecting x and x + i, and there will be a (5 
function in time) current on that link, equal to ji(x) = q a 5{t — t C ross)- If we take the group 
indicies of j to live at the basepoint of the link x, then we should use the value of q which 
departs from site x; the value if we take j to live at the endpoint of the link x + i is the 
adjoint parallel transport of j with group indicies at the basepoint, and hence the charge q 
of the particle when it arrives at site x + i must be the parallel transport of the value at x, 

g(at x + i) = Uj(x)q(at x)Ui(x) . (15) 

If the particle is moving the other direction, from x+i to x, the sign of the current is reversed, 



but (0) still holds. This gives the update rule for q and is the same as was proposed in [30 . 
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Figure 1: An illustration of the degrees of freedom of the proposed system. Lattice sites 
are large dots, and the solid lines joining them are links. Classical fields take values at sites 
(Higgs fields) or on links (connections U and electric fields E). Particles (the small dots) take 
on real valued coordinates and momenta (illustrated with arrows). A particle's charge affects 
the classical fields as if it resided at the nearest lattice site. The dotted lines (really planes, 
extending out of the page) are barriers between the region nearest one site and that nearest 
another, that is, faces of the dual lattice. When a particle crosses a barrier, the charge is 
parallel transported to the new box, the E field on the link orthogonal to the barrier receives 
a kick, and the particle momentum orthogonal to the barrier is changed to conserve energy. 



In what follows we will call such events "boundary crossings" because they correspond to 
a particle crossing the boundary which demarks the volume closest to one lattice site. The 
boundary is a face of the dual lattice, dual to the link in which the current flows. 

We expect that the electric field update rule should be 



dEjjx) 
dt 

Hks 



am 



KS 



dUi(x) 



x.A Z □ 



-TrC/r 



(16) 
(17) 



where i^xs is the Kogut-Susskind Hamiltonian |40| which is standard in the real time field 
literature [H|, and where d/dUi(x) means change with respect to left acting derivatives of 
Ui(x). The appearance of j here is the "meaning" of j and just means that Ei(x) changes 
abruptly by — q a when particle a moves from nearest site x to nearest site x + i, and by q a 
if it goes the other way. 
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Just before crossing 



Just after crossing 





outgoing flux = +q 



outgoing flux = 



outgoing flux = 



outgoing flux = +q 



Figure 2: Boundary crossing, before and after. Before, the flux of E fields (displayed as 
arrows) out of the lefthand box must equal q, and the flux from the righthand box must 
equal 0. After the crossing, the righthand box has flux q and the lefthand box has zero. This 
demands a change of the E field connecting them, but no instantaneous change to the E 
fields which leave these two boxes, since that would change the flux going into other boxes. 



Now let us check whether this definition of E will preserve Gauss' law, that is, whether 
p{x) = Ei{x) -U}{x- i)Ei{x - i)Ui(x - i) . (18) 

i 

(The contribution from the time derivative of U cancels here because it is generated by, and 
hence commutes with, E, and cancels between U and U' .) 
First, we know that 

i oUiix) dUi(x-t) 

which is why Gauss' law is preserved for the classical lattice system without particles. And 
the contribution of j to E cancels p in ( |I8"D precisely because p and j form a conserved 
current. So Gauss' law is indeed conserved. Alternately, the abrupt change to p(x) and 
p(x + i) when a particle passes from being closest to x to being closest to x + i stipulates 
the change to E{(x) so that Gauss' law will still hold at both sites. 

Physically, what this update of E means is that when a particle leaves one box, the 
charge within that box abruptly drops, and the flux of E field out of the box must abruptly 
drop too. Meanwhile, the neighboring box, which received the particle, must have the flux 
of E field out of it abruptly rise. The way to do this is to abruptly change the E field going 
from the first box to the second, by an amount equal to minus the charge, see Fig. |2|. We 
could also do it by changing the E fields along the other links leading out of these boxes, 
but this would spoil Gauss' law at the sites at the other ends of those links and is therefore 
forbidden. Those links will change, in time, via the Hamiltonian evolution of the classical 
Yang-Mills system. 

We should mention that in the abelian theory, one can construct more elaborate Gauss 
law preserving ways to change E, which boil down to moving some of the charge along 
indirect paths between x and x + i. But these are actually not allowed in the nonabelian 
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theory, because the charge transforms nontrivially. To satisfy Gauss' law, the particle charge 
q must also be updated by splitting it up and parallel transporting it along those indirect 
paths; but when they are added together in the final box, they will not all "point in the same 
Lie algebra direction", so the magnitude of the sum will have changed. Similar problems 
arise if one tries to define the charge of a particle as being spread over nearby points, rather 
than residing entirely at the nearest point. Deriving a local, Gauss' law preserving, C, P, 
and T symmetric alternative to our update proposal is highly nontrivial in the nonabelian 
context. 

Before going further we should make a comment about the normalization of fields which 
we have adopted. Since U is a parallel transporter it has been most convenient to use 
the normalization of electric fields in which the Hamiltonian contains a 1/g 2 x E 2 , that 
is our electric field is ga 2 E cont /2. The 1/g 2 in front of the Hamiltonian will be absorbed 
into the temperature, /3l = 4/(g 2 aT), which is customary in the 3-D lattice literature JT7[ . 
Also the particle charge Q 2 will convert into "usual" continuum charge as Q 2 = g A q 2 /4. It 
will be convenient to make the particle momentum appear in the Hamiltonian as = 
Pl(H ks + J2 a Ma), so that it can be directly compared with electric field energies. It is 
then related to physical units by pi att = p C ont x g 2 a/4. Though these normalizations seem 
strange, they scale out all dimensionful quantities and make the numerical theory have the 
least awkward inter-relations. 

It remains to define the update of p. This will be almost uniquely specified by the 
requirement that the system be Hamiltonian. Naively one would expect the influence of the 
electric field on p to be 

Pi = q a E«, (20) 

with E{ chosen to be E on the nearest % type link; but this is wrong, as the system energy 
then has a time derivative of piqEi/\p\, which cannot be removed by a corresponding change 
in E because that would spoil Gauss' law. The energy conserving update 



q a E a_ p -_r 1 _± (21) 



n- qE * ~ p *—i 

does not preserve the phase space measure and is also disallowed; in fact, except at the 
instant when a particle crosses a boundary, measure and energy conservation restrict the 
allowed changes to p to a rotation about some adjoint charged vector B, which must be odd 
under C, P, and T; 

p l = f{p 2 )^kP ] B a k q a . (22) 

So a rotation about the magnetic field is allowed, but not required, to preserve energy and 
the phase space measure. 

However, when a particle crosses the face of the dual lattice separating two points x and 
x + i and induces a current ji(x), it changes the energy in the electric field E^x) by 

A(energy) = A{E 2 /2) = -E^x, t cross - 0) ■ q + q 2 /2 = -E^x, t cross + 0) ■ q - q 2 /2 , (23) 

and p must be changed to balance this energy. If only p^ changes, from Pi^ t to Pifi n , then 
energy conservation is 

P 2 m +P 2 ± + A(E 2 /2) = vW, mit + Pi, (24) 
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which is solved by 

PiM = sign(p iiinit )^/(|p| init - A(£ 2 /2)) 2 - p\ (25) 

if both |p|i n it — A(E 2 /2) > and the argument of the square root is > 0; otherwise there is 
no solution, the crossing is energetically forbidden. If the crossing is energetically forbidden, 
we should set pi^ n = — p^hut, no crossing occurs, no current is generated, and the particle 
turns. Otherwise, p^&n is taken from (|25|) above, the current flows, and the particle crosses 
the boundary. The case where the crossing is from x + i to x follows from this case and 
parity symmetry. 

This proposed update for p conserves energy. What about symplectic measure? The 
measure is clearly preserved at all times that no particle crosses a boundary, since the 
evolution of E, U is Hamiltonian, £ Q depends on p a but not £ a , and q and p do not change 
(except for allowed rotations of p). When a particle reflects, the change Pi —>■ —pi also 
preserves the measure. The only nontrivial case is when a particle crosses a boundary. The 
rotation of q a preserves the measure on LSU(2) and the change E — > E — q preserves the 
measure for E since the change is independent of E, which is defined on a vector space. We 
need only check if the particle phase space measure Hi dpid^i is preserved, that is, whether 

dpifin/ dp jMt dpifijdfy ,init 

,init 



Det 



equals 1. Since the change we are discussing is instantaneous, d^&a/dpj- 
determinant reduces to blocks, 



(26) 
and the 



Det 



Det 



i, fin 



Det 



dpi 



fin 



dp h 



init 



(27) 



dpifin/dp jMt dp iM /d£ jM t 

d£ifin/dp jMt <)0 ::u <KjM -A J ^ Y^jX 

and we need to show that these determinants are inverses. 

Without loss of generality we consider the case in which p\ changes, with p\ > 0. We 
handle the dependence of £ nn on £ init first. Changing ^2,init just shifts the place on the wall 
where the particle crosses, but p is the same at all times; so only £ 2 ,fin is changed, by the same 
amount as the initial change. The same holds for £ 3 . However, changing by d£i changes 
the arrival time at the wall by — rfCi/ u i,init- The particle propagates at rather than -u init 
for c^i/v^init longer than without the change, leading to a change in the final positions equal 
to (vf in - v in u)d^i/vi >iB it, phis the change of d£i in £ 1>fln . Hence, 

t>l, fin M, init 
Det (v 2 fin. - V2,iDit)M, 
. (^3,fin - f3,init)M,: 

Now for the dependence of p nn on p^. Since J92,init = P2,fin and similarly for p 3 , only the 
change in the final value of pifi n need be computed. Taking the appropriate derivatives of 
(H) gives 

l.init / |p|init-A(g 2 /2) ' 
.'l.fln V blinit 

Det 


\pU-A(E 2 /2y 



Det 



i,fin 



j, mit 



,mit 
,init 





1 
1 



^l.fin 

ViMt 



(28) 



Det 



dpi, 
dp 



fin 



j, init _ 



HE 2 /2) P2 

Pl.flnblinit 
1 





A(E 2 /2) P3 

Pl,fin|p|mit 



1 



(29) 



Pi, init 
PlM 



\p\i 
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Now observe that |j5| init — A(E 2 /2) = \p\sn an d that pi/\p\ = V\\ so the two determinants are 
inverses and the measure is indeed conserved. The cases where p<i or p 3 change, or where 
the particle moves in the opposite direction, follow from cubic and parity symmetry. 

If we had allowed p in the orthogonal directions to change, it generically would not 
preserve the measure, so the choice of update is almost unique.0 

The first four conditions for a valid update rule have specified the update uniquely except 
for the freedom to perform rotations of the form shown in fl22|). It would seem natural to 
define the magnetic field at a point in terms of the 3 nearest plaquettesf] and to rotate 
momenta according to 

iH = e ijk ?i-Btf. (30) 

As we discussed, the magnetic term in Wong's equations plays no role in reproducing the 
right hard thermal loops for a plasma close to equilibrium. Thus we are free to omit this 
rotation when dealing with quasi-equilibrium processes. 

Next let us verify that the continuum limit of our update rules give Wong's equations. 
This is known for the classical field equations without particles evolving under the Kogut- 
Susskind Hamiltonian, so we need only check the terms involving particles. The update 
of the particle charge is explicitly the adjoint parallel transport along its trajectory, as it 
should be. The current in the % direction from a particle occurs in jolts of magnitude q and 
frequency 1/vi, so the time averaged current is qvi as it should be. The impulse on a moving 
particle is of magnitude Api = p i fin — Pi^t, expanding ( |23"D and fl25|) to leading order in q 
gives Api = q-Ei/vi, and such impulses also occur with frequency so the time averaged 
force on the particle is q ■ Ei, also as it should be. We will explore the corrections to the 
plasma frequency, due to the discrete nature of the current and the "kicks," in Appendix A. 
Wong's equations are recovered for the motion of a particle through slowly varying fields, if 
q is sufficiently small. 

We still must check whether the thermodynamics of this system are correct; we will do 
this in the next chapter. But first, let us review and summarize the update rule. The fields 
evolve under the Kogut-Susskind Hamiltonian and the particles move freely, except at such 
exact instants when a particle crosses a boundary (a face of the dual lattice), i.e. when it goes 
from being nearest one point to nearest another. The charge is then parallel transported by 
the link operator dual to the face, and the electric field on that link is abruptly changed by 
— q (or +q if the motion is in the — direction). The particle momentum parallel to the link is 
changed to cancel the energy change of the electric field. If that is energetically impossible, 
the particle momentum parallel to the link is flipped and no crossing, or change to q, E, 
occurs. This update is Hamiltonian with total energy 

Energy = £ + Wl - W) + £ M« , (31) 

x,i A □ V Z J a 

and preserves Gauss' law identically. It is also manifestly gauge, C, P, T, and cubic invariant. 

4 A rotation of p± violates parity unless its magnitude depends on B or some other P odd field; and an 
energy conserving change in the magnitudes of p\ and p± will have nonunity Jacobian. 

5 One should be careful here to make the choice of plaquettes in a way which is preserved under parity 
and the cubic point group. 
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4 Thermodynamics of the lattice system 



In the last section we have proposed a continuum time system and its update rule and 
have shown that it is Hamiltonian, ie that the update rule preserves energy and phase space 
measure. It therefore has well defined thermodynamics, which we now explore. The canonical 
partition function is 

Z = f]ld 3 ^ a d 3 p a dq Q ]ldE i (x)dU i (x)]l6(-p(x) + D L ■ E(x))e~^ H , (32) 

a i,x x 

h = E^+Efi-W+Ebla- (33) 

x,i Z □ V Z J a 

The appropriate measures, ie the Haar measure for U and the Lebesgue measure on LSU(2) 
restricted to q 2 = Q 2 for q, are implied. The 5 function enforces Gauss' law at each point 
and the meanings of p and ■ E are as previously defined. The quantity /3l here combines 
all the dimensionful parameters of the lattice system, /3l = 4/g 2 aT; when it is large we are 
on a fine lattice, or at weak coupling, which is equivalent in the classical theory. 

Note first that the kinetic momenta will obey a Boltzmann distribution and are indepen- 
dent of all other degrees of freedom; the partition function factorizes and if we are interested 
in the thermodynamics of the IR classical gauge fields then we may integrate out p. Note 
also that £ and q do not appear in the Hamiltonian; their only influence is in determining p 
appearing in Gauss' law. 

It is most convenient to enforce Gauss' law with an adjoint valued Lagrange multiplier 



A [17 



Y[(D L ■ E(x) - p(x)) = f J] dAo(x) exp 



i(3 L J2A a (x)(D L -E a (x)-p a (x)) 



(34) 



which makes the electric fields quadratic; the E integral may now be performed, generating 
a kinetic term for Aq. The partition function reduces to 



UA 



HdU i (x)dA (x)e-^ HuA I(A ) , 

i,x 

Wl-iw D ) + i;OAA) 2 (*) 



I(Aq) = J l[d 3 £ a dq a e- ll3 ^ Ao{x)pix) 



(35) 

(36) 
(37) 



where the meaning of (D-^Aq) 2 should be clear. The particles enter in the last term, which 
only depends on A . 

Now let us compute I(A ). Since p is the sum of p a from each particle, exp — i[3 L A p 
is the product over a of exp —i(3 L A p a ] the integral factorizes into an integral over each 
particle, 



i(a ) = n 



J dq a J d 3 £ a exp (-i!3^A a (x)q a a 



x 



1 £ at x 
£ not at x 



(38) 
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The integral over £, normalized so f d£ = 1, gives a sum over sites of a term where the 
particle is at that site, 



km = n 



Now using 



(igexpH/3 L y^(a;)g a ) 



sin (^Q^ 



(normalizing so that / dq = 1) and performing the product, we arrive at 

" 1 sin (> L Q V / A^1 

y ( A oJ- 2^ 



(39) 



(40) 



AT3 



(41) 



The above expression is exact but not very insightful as written. It is best to take a 
thermodynamic limit and to expand the expression in the sum. Denote by (ARG) the mean 
value of the argument of the sum in (f4T|). Rescale I{Aq) by (ARG) - ^, which just changes 
the normalization of the partition function, and write it as 



I (A 



1+ -1 + 



E 



sin (faQyfAl 



A 3 (ARG) x PlQ JaJ^ 




(42) 



As A 3 — > oo and N p — > oo with N p /N 3 = (n) fixed, the term in round parenthesis vanishes 
as N~ 3 / 2 , and we may use the approximation 



1 + 



x 
N 



N 



exp(x)exp(-x 2 /2N) x (1 + 0(x 3 /N 2 )) 



(43) 



The exp(— x 2 /2N) term in the identity means that there is a very weak nonlocal interaction 
term between fluctuations from the mean value of Aq at pairs of points, trying to force the 
global average of Aq towards its equilibrium value. It is probably safe to ignore this term, 
and higher corrections strictly vanish in the large N limit. Neglecting the nonlocal term, we 
get 



I(Aq) ~ exp 



(n) sin (PlQ/a 2 



(ARG) 



PlQJa 2 ( 




(44) 



Since we are typically interested in a system where A 3 ~ 10 4 and N p > 10 5 , the thermody- 
namic limit is justified. 

Now, expanding the term in the sum, 



sin lp L QJA 2 (x 



1 PlQ 2 A 2 (x) | PtQ\A 2 ) 2 (x) 



PlQJA&x) 



6 



120 



(45) 
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we see that I{Aq) contains a Debye screening term and quartic and higher self-interaction 
terms for A ; I(A ) becomes exp(— P l V(Aq)), where the potential V is 



< = (47) 



3(ARG) 
(n)(3jQ A 
30(ARG) 



(48) 



Finally, 



(ARG) = 1 - @f-(A*) + O(PtQ 4 A* ) . (49) 

(Aq) is UV dominated and well approximated by its perturbative value. At lowest order, 
(Aq) = 3S/47t/3l, with £ ~ 3.1759 |p!0| . In practice it is necessary to make Q ~ P^ 1 , see 
below, in which case, for /3l ~ 10, (ARG) typically differs from 1 by ~ 1%. As we men- 
tioned in Section |, the ratio m^/ (n) should equal (2qf)(E^ 1 ) = (Q 2 /3)(2/5l/2) for classical 
particles in the continuum. Discretization has shifted this by a small 0(Q 2 ) correction. 
In a group other than SU(2), the calculation begins to go differently at Eq. (f40"D, and 



the final Debye mass is different. For instance, in U(l) gauge theory, the integral there gives 
cos(/3l£*Mo)> and m\> and Aa are 3 and 5 times larger, respectively. 

We have now shown that the thermodynamics are the same as lattice regulated Yang- 
Mills theory in the dimensional reduction approximation, at a specific value of Debye mass, 
except for the added higher order interaction terms for the Aq field. A minimal requirement 
for the thermodynamics to be all right is that the Aq field is not strongly coupled, Aa < 1. 
The lattice Aa corresponds to the physical AXx/g 2 , so this requirement is that the coupling 
of the Aq sector is weaker than the coupling of the gauge sector, which should be sufficient 
since the Aq field is quite massive. This requirement is roughly 

30 

Q\n)<^. (50) 



Pi 



In fact it will turn out that dynamic considerations require that Q be on order or smaller 
than in which case, unless (n) is very large, Aa will be very small, as it should be in 
the dimensionally reduced Hamiltonian. 

We should also require that the Debye screening mass from particles be larger than that 
from hard lattice modes (as otherwise the hard thermal loops are dominated by wrong lattice 
mode contributions, rather than right particle contributions), in which case 

1 1Y, 

-M<(n)Q 2 . (51) 

However, it is not necessary on thermodynamic grounds to make small in lattice units. 
When m\y is large in physical units, ie m\ 3> (3^ 2 * n lattice units, then the influence of the A 
field on the vector fields is perturbative. Making m\, on order the lattice spacing complicates 
the problem of integrating out the Aq field to find its influence on the infrared physics, but the 
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one loop integration can still be performed; we treat this problem in Appendix B. However, 
it will turn out on dynamical grounds that one should not make the Debye mass too large, 
as the plasma frequency then takes on lattice artifact corrections, discussed in Appendix A. 

5 Discrete time update algorithm 

In this section we give a detailed description of a stable and accurate discrete time update 
algorithm to solve the equations of motion for the particles and lattice fields. Probably most 
readers can safely skip this section; nothing in it is essential for understanding the rest of 
the paper, although it certainly is important that a stable and well behaved algorithm exists 
for implementing the system described in the last two sections. 

The numerical algorithm we will construct is essentially a leapfrog, with modifications 
to meet the needs of the specific problem under investigation. In particular we make sure 
the algorithm is time centered, and the particle update part is exactly energy conserving. 
The stepsize errors in the algorithm should be 0((At) 2 ) where At is the time step in lattice 
units. 

The variables to be updated are: (£ a , p a , q a ) and (E x>i ,U x>i ). If the connections U were 
fixed and the particles did not rotate due to magnetic fields, then we could perform the 
update (of E, p, q, £) exactly, as follows. Starting with E, p, £, q at time t, we project where 
£ will be at time t + At, 

£proj,a,i(t + At) = ^i(t) + Atff- . (52) 

\P\a 

For all a where the projected location is in the same box (dual lattice cell) as the initial 
location, ie no boundary is crossed, then £(t + At) equals the projection, and we update 
£. For other points, we draw the straight line path between £(£) and £ pro j and find the 
first boundary it crosses. £ can be updated to this point, Abound) which it reaches at time 
Abound = t+ Abound — At this point we must be careful, because more than one particle 
may interact with the same electric field, and the updates will depend on the order of 
interaction. So we order all particles which cross a boundary according to crossing time; 
then we solve the crossing conditions in that order, modifying E, q, and p as discussed in 
Section |3|. For each particle, after these variables are updated, the particle's position at 
time t + At is again projected, starting at Abound and using the new p. If it crosses another 
boundary before time t + At, we again update it to the boundary, compute the time of 
arrival, and insert it at the place appropriate for the new crossing time in the collection of 
particles to be updated; otherwise we update it to its projected position at time t + At. 
When the last crossing has been dealt with, then all E, p, q, and £ have been updated to 
time t + At. This algorithm is exact except for roundoff errors. 

One comment is in order about this update. Ordering the particles by time of crossing 
takes 0(N p \n N p ) steps, so the algorithm does not quite scale linearly with volume. In 
practice, though, the ordering takes CAtN p In N p computations with C a fairly small number, 
and this part of the algorithm takes less time than formally 0(N P ) or 0(N 3 ) parts. 

Now we must incorporate the above idea into the leapfrog update of the E, U fields. The 
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leapfrog in the absence of particles is 



Ui(x,t + At/2) = exp(iAtT a E?(x,t))Ui(x,t-At/2), (53) 
^, t + A<> - *M)-* ^' + + *g» . 0-) 

where the meaning of OH-^/dU is defined between ( ]T6| ) and QTST) . If one needed to define 
Ef{x, t + At/2) one could do so by applying only half of the update, fl54|) . 

We combine this leapfrog and the update of £ etc. discussed above, as follows. 

1. Start with E(t), p(t), q(t), f (t), and U(t - At/2). 

2. Determine U(t + At/2) according to (|53|). 

3. Apply half the i£, Z7 leapfrog update of E; namely, set 

„/ s AtdH KS (U(t + At/2)) , . 

Also update to p(t + 0) by rotating about B for time At/ '2; first set 

Ptemp.a.t = Pa,i + ^ e ijfcf^p^fc (0<f . (56) 

(Here B k (£) means the i,j plaquette which is closest to the point £ and with indi- 
cies living at the point where the indicies of q reside.) Then rescale p to its original 
magnitude, 

Pa,i(t + 0) = Ptemp,a,i , ^ . (57) 

IP | temp, a 

4. Update E, q, £, and p from time t + to time t + At — using the "fixed connection, 
no B field" algorithm presented above, and U = U(t + At/2). 

5. Apply the other half of the E, U leapfrog, 

and rotate the momenta as in step 3. Now, return to 1., but with the value of t 
incremented by At. 

Applying these in order is one leapfrog update. Note that the update is time symmetric 
and exactly Gauss constraint preserving at each step, and that step 4. exactly conserves 
energy. The overall conservation of energy is exactly as good as in the Kogut-Susskind 
leapfrog algorithm, that is, energy fluctuates by a small 0((At) 2 A^ 3 / 2 ) amount (which for a 
20 3 grid at At = 0.05 and 30 particles per site is less than a part in 10 5 of the total energy) 
and the central value of the energy is absolutely stable. 

As we have discussed, the effect of magnetic fields on particles does not contribute to the 
hard thermal loops for a plasma close to equilibrium. So that part of steps 3. and 5. can be 
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left out when studying quasi-equilibrium processes, which in practice saves at least 1/3 of 
the update time. 

Finally we present a canonical ensemble thermalization algorithm for this system. It is 
a straightforward generalization of the "constrained molecular dynamics" algorithm of |[42| . 
Beginning from an arbitrary choice of [/(At/2), £(0), and q(0), we choose p(0 + 0) from the 
Boltzmann distribution at some inverse temperature /3l and E(0 + 0) from the Gaussian 
distribution at the same temperature, modulo the Gauss constraints. As in [52] this is done 
by choosing E without regard to the constraints and then orthogonally projecting to the 
constraint surface, which correctly thermalizes the transverse components of E and correctly 
enforces the constraints on the longitudinal components. The algorithm is identical to that 
42| except that the particle contribution to Gauss' law must be added. This chooses 



m 



p, E with thermal weight from the fixed U, £, q subspace of phase space. Then we evolve 
the system under the Hamiltonian evolution, using the algorithm presented above, for some 
length of time, at the end of which we again draw E and p from the thermal ensemble; we 
repeat until the (athermal) original information in the U, £, and q has been destroyed and 
measurables attain values which do not change in the mean under Hamiltonian evolution or 
further thermalization. The spirit of this algorithm is that of a molecular dynamics Monte- 
Carlo. Note that it is important to choose E and p at some time between t + and t + At — 
and not at t or t + At because the half update from t — to t changes E from being Gaussian 
and uncorrelated with dH^/dU to being correlated with BH^/dU . (Similarly in |42| it was 
necessary to choose E defined at the same time as U and to perform a half update before 
beginning the leapfrog.) 

We should also note that there is no obstacle to applying a Langevin type thermalization 
algorithm based on the one developed in and in particular that it is trivial to couple 
Langevin noise to the particle momenta. This might be important if one wanted to simulate 
thermalization of these modes through interactions with some other degrees of freedom, for 
instance strong scattering of fermions in the electroweak model. 



6 Some results for the abelian theory 

Before diving into the study of Nqs diffusion we should check that the particle method is 
producing the right physics of HTL's. In this section we will first discuss how small Q 2 
must be for the system to give good behavior, and then we will study the abelian theory at 
suitably small Q 2 to see if the hard contributions to the retarded self-energy are correct. 

We saw in Section |] that a necessary condition for the theory to have a weakly coupled 
A sector is 

3D 

Q»<-3. (59) 
Pl 

In fact, dynamical considerations demand that Q 2 be still smaller. To get the right hard 
thermal loop effects in the dynamics, we need a particle to travel a distance longer than 
the magnetic length scale l/g 2 T (or /3l in lattice units) before its momentum is randomized 
by interactions with the plasma. Normally we would expect the dominant randomizing 
process to be Coulomb scattering from other particles as depicted on the left-hand side of 
Figure |3|; but in fact the dominant processes are absorption or emission of hard classical field 
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Figure 3: Diagram which naively dominates particle scattering, left, and Cerenkov diagrams 
which are really dominant, right. 



excitations. These processes are not kinematically forbidden because the lattice dispersion 
relations "turn over" at high momentum, so hard excitations move slower than the speed 
of light, allowing particles to Cerenkov radiate or absorb. From the point of view of the 
lattice modes, this is saying that the ultraviolet lattice modes are Landau damped due to 
the particles, which is possible because their dispersion relation has u/k < 1. 

A "worst case" estimate is that the typical electric field which a particle interacts with 
is uncorrelated with the charge q of the particle, so the momentum p falls by q 2 /2 = Q 2 /2 
in addition to receiving a kick uncorrelated with its initial value. The average particle 
(averaging over directions of motion) crosses walls with frequency 3/2 in lattice units; so the 
decorrelation rate for p under the above approximation is 

rate ~ W) = — ' (60) 

which must be <C P^ 1 ; so Q 2 <C 4//3f\ In practice the above estimate for randomization 
is close to the real behaivor; for /3l = 10 and Q = 0.08, the decoherence time is about 80 
in lattice units. In our work we typically demand Q 2 < j3^ 2 , which forces the number of 
particles to be quite large, (n) ~ 30, in order to satisfy (0). This makes the update of 
particles the dominant numerical cost, but it pushes us closer to the Vlasov equation limit. 

Note also that the processes mentioned above make the evolution of the UV classical 
lattice modes damped and noisy, with a damping strength proportional to Q 2 (n). The 
consequences for the infrared magnetic sector deserve investigation. 

The above behavior is disturbing enough to encourage us to check that the behavior of 
the infrared degrees of freedom is correct. This is most easily done for the abelian theory, 
because there the electric field and the current are gauge invariant and one can easily probe 
the system with external currents and study the response. It would also be possible to do 
this in the nonabelian theory in a specific gauge, but complications from pure classical gauge 
theory interactions complicate interpreting results; the abelian theory provides a nice, clean 
environment to study the dynamics of the particle technique. 

We will perform two numerical tests on the abelian theory. First, we verify that the lattice 
field has the correct dispersion in the infrared. Second, we probe the retarded propagator 
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Figure 4: Dispersion relation: free field in the continuum (dashed line: to — k) ; free field on 
the lattice (dotted curve: ua = 2 sin (/ca/2)); plasma modes: transverse (theory: the upper 
solid curve; data: solid rectangles) and longitudinal (theory: the lower solid curve; data: 
open rectangles). 

by studying the linear response of the abelian plasma to an external current. 

To measure the plasma dispersion relation, we pick initial conditions so the (spatial) 
Fourier mode of wave vector k has a large electric field. This sets up a plasma oscillation 
with this wave number. We evolve the fields and measure the electric field in the same 
Fourier mode at regular time intervals. Then we Fourier transform this time series and pick 
out the frequency with the most power; this is the plasma frequency. We have used this 
procedure to measure the dispersion for both transverse and longitudinal plasma modes. 
The results are summarized in Figure (|, where we compare them to the lattice, free field 
dispersion relation and the continuum dispersion relations, with and without hard thermal 
loops. These data are for Q = 0.0159, /3l = 18.86, and (n) = 50, corresponding to a plasma 
frequency of u p = 0.282. (Recall that in the abelian theory and uo 2 are 3 times larger in 
terms of Q 2 (n) than in SU(2).) The numercal results agree with theory remarkably well in 
the infrared. As k gets bigger, the data points deviate from the continuum theory curve and 
bend down to match the lattice dispersion relation. This is expected because hard plasma 
modes are much less influenced by hard thermal loops and behave like free lattice modes. 
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We have also computed the plasma frequency at k = analytically, in Appendix |A|. There 
we conclude that the correction due to lattice artifacts should be negligible for u£ <C 1 in 
lattice units. This is supported by our numerical results. 

To probe the retarded (transverse) photon propagator, we drive the system with an 
external current of the following form: 

ji = jijo sin (ut) sin (kx) , (61) 

where the amplitude jo is small so that linear response theory applies. We study response 
of the plasma and measure the space-time average of j ■ E, which, according to Appendix [C], 
can be written in terms of the transverse photon polarization function: 

/ • E) = Jo uUj(uj,k) 

U ; 4 [u 2 -k^ -Tl r {uo,k)f + Tl 2 {uo,k)' 1 ' 

where H T (u,k), Hi(u,k) are the real and imaginary parts of the transverse polarization 
function, respectively. We also measure the out of phase response of E and the case of 
longitudinal excitation; the specific expressions are in Appendix |C[ 

We want both to compare to theory and to test the rotational invariance of the particle 
induced hard thermal loops. To do so, we choose two k vectors which are of the same length 
but are inequivalent under the cubic point group, and we study the complex propagator 
for each, at frequencies above and below the plasma resonance. We choose the vectors 
k = (3,0,0) x tt/16 and k = (2,2,1) x vr/16 on a 32 3 lattice, with f3 L = 20, Q 2 = 0.0005, 
and (n) = 60 (so u p = 0.447). We excite the plasma at each wave vector and numerous 
frequencies in turn and integrate the resulting response for long enough to get reasonably 
clean results. 

We present the results, plotted against the theory, in Figure [5]. A few comments are in 
order. In the strict hard thermal loop approximation there would be no imaginary part to 
the self-energy above the light cone, but the appearance of a small imaginary part, which will 
arise at 0(Q A (3^(n)) since we are in the abelian theory, qualitatively changes the response 
near the resonance by giving the resonance some width. For the transverse data, we have 
added a small phenomenological imaginary part to the self-energy chosen to make the theory 
give a resonance of about the same sharpness as the data. The match between data and 
theory is very good except that the location of the resonance is shifted a little towards lower 
uj in the data. This is expected; the particle response weakens with frequency because the 
particles cannot respond to a field faster than they encounter it, which is how often they 
cross a boundary, as we examine quantitatively in Appendix A. Also note that the data for 
k oc (2,2,1) (triangles) have their resonance very slightly after the data for k oc (3,0,0), 
which is expected from the 0(k 4 ) lattice corrections to the dispersion relations. However, 
besides this, they are indistinguishable, which is a good check that the particle contribution 
to the self-energy is rotationally invariant. 

We have not added a small imaginary part to the self-energy above the light cone in the 
theory lines in the plots for the longitudinal propagator, but the data clearly show that one 
is present. The agreement between data and theory is generally good for both longitudinal 
and transverse propagators, and the data for the two values of k disagree by about the same 
amount as the jitter in the data caused by statistical error. The comparison to the theory 
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Figure 5: Retarded Propagator-actually iuD R . Above left, the out of phase (nondissipative) 
part of the transverse propagator; lower left, the in phase (dissipative) part of the transverse 
propagator. The right figures are the same for the longitudinal propagator. In each case, 
the solid lines are theory, the squares are the k = (3,0,0) x n/16 data, and the triangles 
are the k = (2, 1, 1) x n/16 data. In the transverse theory figures, the curves with smaller, 
filled data points are rescalings of the data and theory, so the resonance will fit in the plot; 
the rescalings are by 8 and 120 for the out of phase and in phase responses respectively. All 
numbers are in lattice units. 
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without hard thermal loops is stark; for instance, the in phase response would then be zero 
in all cases, and the out of phase, longitudinal response would have no resonance but would 
behave simply as 

We conclude that the particles add hard thermal loops which, for k <C ir in lattice units, 
are rotationally invariant and very close to the correct "hard thermal loops", except for 
corrections which are O(uo 2 ) and 0(Q A (3^(n)) (both in lattice units). This is totally unlike 
the hard thermal loops induced by hard lattice modes in the nonabelian theory or the abelian 
Higgs model, which are rotationally non- invariant [p0|. 



7 Diffusion of Chern-Simons number 

An outstanding question which the method developed here can answer is; what is the diffu- 
sion constant for Chern-Simons number Acs in the symmetric electroweak phase, or in pure 
Yang-Mills theory? In the continuum, the Chern-Simons number is defined as 

2 fiva(3 

PK, = ^—1*1%, (63) 
A cs = / d 3 xK = JL; J d 3 xe l3k (f«AI - 9 -f abc A<tA)A^ . (64) 

where Latin indices (ijk) imply sums over the three space directions with positive metric. In 
vacuum Acs is an integer, equal to the winding number of the gauge transformation which 
carries the gauge fields into the trivial fields A = 0. This can be nonzero for groups with 
nontrivial third homotopy group. Chern-Simons number is of interest physically because of 
how it is related to the behavior of fermions coupled to the gauge fields; a vacuum to vacuum 
process which changes Acs by n pulls n left handed negative energy solutions of the Dirac 
operator up to positive energy and pushes n right handed solutions from positive to negative 
energy, if the fermion couples to the gauge fields in the fundamental representation. In a 
chiral theory, like the SU(2) sector of the standard model, no right handed particles couple, 
and there will be net particle creation (recall that in vacuum, negative energy solutions are 
occupied and positive energy solutions are not). Summing over the A F = 3 generations one 
finds that the baryon number and lepton number currents both have nonzero divergences, 

d% L = d% B = A F «^A M , (65) 

so the number of baryons Ab and of leptons Al changes as 

N B (t) - A B (0) = N L (t) - N L (0) = N F (N cs (t) - N cs (0)) . (66) 

At high temperatures the efficiency with which Nb is violated is related to the diffusion 
constant for Acs, 

r.iimiim (WsM-iW) (67) 
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(where () refers to the thermal ensemble) by a fluctuation-dissipation relation fH and 



standard thermodynamic arguments ||. In the minimal standard modelQ, one finds 

1 d(N B + N L ) 39 



N B + N L dt ~~ 4T 3 



r . (68) 



It would be phenomenologically interesting to know T at high temperatures, which pre- 
sumably existed in the early universe before the electroweak SU(2)xU(l) symmetry was 
spontaneously broken. 



There has been a great deal of work to date on determining T |L7|, [15|, [41], ^2j| , but 
recently there have been two important developments. 

The first is an analytic argument due to Arnold, Son, and Yaffe (ASY) [p3fl . They point 
out that a change in Nqs involves the evolution of very infrared magnetic fields. The hard 
thermal loops cause these fields to evolve in an overdamped manner. This is familiar from 
the study of abelian (electromagnetic) plasmas; magnetic fields of wavelength A ^> l/w p 
get "frozen" by the conductivity of the plasma and evolve on the time scale r ~ A(Au; p ) 2 , 
assuming A is much shorter than the diffusion length of the charge carriers. ASY argue that 
the same physics should apply in the nonabelian plasma, at least for A less than or on order 
1 / g 2 T; so T should be parametrically of order 

^M^) (69) 

If this argument is correct then T depends strongly on the physics of hard thermal loops, 
which we must get right to find the correct T. Classical lattice theory by itself (without 
classical particles) does not. The damping coefficient describing the overdamped evolution 
of infrared magnetic fields on long time scales grows linearly with 1/a and is not rotationally 
invariant [55]. 



The second development is that better definitions of Nqs on the lattice have been devel- 
oped J35J, [25L 5B|. Previous definitions in terms of local operators contained lattice spacing 



dependent systematic errors. Moore and Turok proposed a definition which is topolog- 
ical and hence avoids such errors. Their results verify that T depends strongly on lattice 
spacing, in a manner which appears consistent with the ASY scaling law, though only if 
there are substantial corrections to that law at higher order in (g 2 T) 2 /uo 2 . 

By using the topological method of Moore and Turok to track the evolution of Nqs, 
and by using the classical field theory with particles added to correctly reproduce the hard 
thermal loops, we can now get a determination of the diffusion constant for Ncs which 
accounts correctly both for topology and for hard thermal loops. We can also ensure that 
the thermodynamics of the system under study are correct by using the 0(a) improved 



matching developed in |Tl[ and extended to arbitrary Debye mass in Appendix ||. 

We will do so in pure Yang-Mills theory, which should correspond to the very high tem- 
perature limit of the Standard Model because the thermal Higgs mass becomes large enough 
at high temperature that the Higgs field can be integrated out (though its contribution to 

6 In extensions with light baryon number or lepton number carrying particles, 39/4 will be replaced by 
something smaller; in the supersymmetric theory, if all the squarks and sleptons were light, 39/4 would 
become 39/12. 
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the hard thermal loops should of course be included). The sphaleron rate in the symmetric 
phase at the equilibrium point of the phase transition will differ somewhat from the Yang- 
Mills theory value, in a way which depends on the (unknown) couplings of the Higgs sector. 
It is straightforward to add the Higgs field to the theory we have developed-in particular the 
influence of hard modes on the Higgs field can be completely accounted for by the choice of 



Higgs mass |47[-so there is no obstacle to extending what we do here to that case. 

There is an important complication to our plan; in a nonabelian theory, the UV classical 
field modes will also generate hard thermal loops; and as Bodeker et. al. |2(J have shown, 
the functional form of the hard thermal loops they provide is NOT the same as Eq. (|TT|). 
Hence, the actual hard thermal loop contribution to the lattice system with particles will be 
of form 

m D,particles ' particles [A] + ^D,UV lattice ' ^UV lattice [A] , (70) 

where r partic i es [A] has the correct and r uvlattice L4] has the wrong functional form. This is 
further complicated because of the interactions between the particles and the UV lattice 
modes. The UV lattice modes are Landau damped, and as we increase m 2 ^ from particles, 
that damping becomes stronger. If this damping is strong enough, then UV classical modes 
have short propagation distances, and they will not propagate interactions over large spa- 
tiotemporal separations. The hard thermal loop effects most important to the iVcs diffusion 
rate are those on length scales of order l/g 2 T ||23|| , so the relevant contribution from UV 
classical lattice modes may have an m 2 ^ dependent suppression. We do not know a good way 
to estimate the importance of this suppression, so we will treat it as a source of systematic 
error. 

We want to know the result of a triple limit. The innermost limit is the limit of Q 2 — ► 
and (n) — ► oo with Q 2 (n) fixed; in this limit the particles generate only hard thermal loop 
effects. The next limit is the limit as ^d, particles / m D,uv lattice becomes large; in this limit 
the hard thermal loops are of the correct functional form, plus a correction which is small 
relative to the total strength. At fixed lattice spacing, this limit means making large; 
so we can only learn about the parametrically leading behavior in the limit of large mp. 
Finally, we should take a small lattice spacing limit. 

Limited numerical resources make it impossible to really achieve this triple limit. We 
will choose a value for Q 2 small enough that we can expect to be in the relevant limit there, 
and we will also assume that (3l ~ 10 puts us far enough in the small a limit, if 0(a) 
thermodynamic corrections are usecf]. We will check both of these limits by varying Q 2 
holding Q 2 (n) fixed and by varying a, to verify that the dependence is weak; but we make 
no serious attempt to extrapolate to these limits. We concentrate on what we consider the 
most phenomenologically interesting limit, what happens as we make m 2 ^ large. In particular 
we want to know whether the Yd scales according to the ASY prediction, Yd oc m^ 2 . 

We will then try to check three things: 

1. T, in physical units, should depend weakly on the lattice spacing, provided that the 
physical value of the plasma frequency is held fixed; 



7 We verify in Appendix [a| that the frequency corrections to the particle hard thermal loops are 0(u 2 a 2 ), 
and we believe the same should be true of the finite k corrections. This leaves nonrenormalizable operators, 
also 0(k 2 a 2 ), thermodynamic errors beyond 0(a), and an O(a) rescaling of the time scale, which we estimate 
in Appendix |5] but have not computed. 
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Figure 6: Nqs during a Hamiltonian evolution, tracked by the slave field method (which 
returns the integer value of Ncs for the nearest vacuum). This is part of the data for row (5) 
of the table. The back and forth stuttering accompanying many winding number changes is 



an expected feature if the ASY picture is right [23 1 



2. T should depend on Q 2 and (n) only through the combination Q 2 (n), when we have 
chosen Q 2 small enough that the particle trajectories are ballistic on the nonperturba- 
tive scale (g 2 T)~ l ; 

3. T should depend inversely on hard thermal loop strength, up to corrections due to the 
UV classical lattice mode contributions. 

To test 1., we measure T at three lattice spacings, but with the same value of Q 2 and (n) 
in physical units. To test 2., we double the value of (n) and halve the value of Q 2 relative 
to the runs used in 1., and to test 3. we double and quadruple (n) relative to the runs used 
to study 1., but keep the same value of Q 2 . If the ASY scaling law is right, then up to 
corrections due to HTL's from UV classical lattice modes, all these runs should give the 
same value for T, except for the last, where T should be about half and then a quarter as 
large. If T does not depend on hard thermal loops, then all the results should be the same, 
and should be the same as the result without any particles added. 

Our results are presented in Table [I]. Each data point is extracted from several Hamil- 
tonian evolutions from independent thermal initial conditions. We present a sample of iVcs 
during such a Hamiltonian trajectory in Figure ||. The analysis techniques used to extract 
T are the same as in For each choice of parameters the sum of lengths of evolutions 
is about 90000 lattice units of time. The value of (n) for the /3l = 8, 10, 12 data keeps the 
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Table 1: N cs diffusion constant in physical units, Y = Ka^T 41 , varying lattice spacing, 
particle charge, and particle number. The inverse lattice spacing /3l = 4/g 2 aT used here 
is the one including the perturbative corrections found in [|TT| and Appendix ||. The last 
column is the coefficient of the ASY scaling law, see the text. 



physical particle density, which is proportional to (3^ (n) , fixed; and as we have discussed, Q 2 
does not scale with lattice spacing. Hence these three results test (1). Row 4 has the same 
value of Q 2 (n) as row 2, so comparing them tests (2). Finally, comparing rows 2, 5, and 6 
tests (3). 

The results are expressed through the dimensionless quantity k, which in the continuum 
is defined through 

T = Ka 4 T 4 , (71) 

or on the lattice, 

T = k((3 LMp 7t)- 4 (72) 

times a correction, discussed in Appendix 0, to account for the correct matching of the time 
scales. 

We also present the results in terms of the coefficient of the ASY scaling law, which we 
write as 



Y = k' 



a 5 T 4 . 



mt 



(73) 



We use here because it most conveniently characterizes the size of hard thermal loop 
effects. For the particle degrees of freedom, the Debye mass in physical units is 



g 2 (n)(latt. units)^ 4 T 2 
48 



(74) 



which for the first 4 columns equals Ag A T 2 , a little less than the physical value, which is 
ll 9 2 T 2 /6 )9 2 ^0.4 |. 

Really the ASY scaling law says the results should depend not on the Debye mass but 
on a damping coefficient proportional to the transverse self-energy at u ^ k ~ g 2 T. This is 
simply related to m 2 ^ in the case that the hard particles have a rotationally invariant spectrum 
and move at the speed of light. The particle degrees of freedom satisfy this requirement, 
but the hard classical lattice modes do not, so we have accounted for their contribution 
using the techniques of Arnold P5fl . His result is that the ratio of damping coefficient to 
Debye mass squared is roughly (0.68 + 0.2) times smaller for hard classical lattice modes than 
for ultrarelativistic particles. However, Landau damping of the UV classical lattice modes, 
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mentioned earlier, may suppress their contribution to the transverse self-energy at k ~ g 2 T; 
since we do not know how to compute the extent of this suppression, there is a systematic 
error. The upper limit of the systematic error bar we present is if they contribute fully, in 
which case we add 0.68 times the Debye mass squared from classical lattice modes to that 
from particles when converting k to k'. (In physical units the Debye mass squared from 
classical lattice modes is mj^latt) = (T l f3j J /4:7r)g 4l T 2 .) The lower limit is if their contribution 
is fully frustrated by Landau damping off of hard particles, in which case we just use m 2 ^ 
from particles to convert from k to n'. The larger Q 2 (n), the stronger the Landau damping; 
so the systematic is not common to all runs. 

Our results are roughly consistent with lattice spacing independence and with dependence 
on Q 2 and (n) only through the combination Q 2 (n) . There does seem to be a weak systematic 
trend in lattice spacing, which could be partly from 0(a 2 ) effects we have not attempted to 
compute. For instance, besides the 0(a) corrections to the thermodynamics computed in 
|TT|| and Appendix [FJ, there are two 0(a 2 ) corrections: a renormalization of the coupling, and 



a nonrenormalizable (-Dj-FV,) 2 term, which appears in the Hamiltonian with a negative sign. 
Both would raise the rate on coarser lattices. There are also 0(a 2 ) corrections to the Wong's 
equation limit of the interactions between particles and long wavelength modes. There also 
may be a weak trend in Q 2 when Q 2 (n) is held constant, because we are not sufficiently 
close to the small Q 2 and large (n) limit. Not being in this limit means that the eikonal 
approximation used to turn the particles into hard thermal loop effects is not quite true. 
Scattering of the particles will tend to reduce their effectiveness. Hence one might expect a 
weak systematic where k rises with Q 2 at Q 2 (n) fixed. Row 4 in the table suggests this but 
the effect is not very statistically significant. Systematics have not been eliminated but they 
are small. 

Our results rule out the "old picture" that T should depend on a 4 T 4 by demonstrating 
the importance of hard thermal loop effects. If the a 4 T A law were correct, all the values for 
k would agree, and would agree with the value in Yang-Mills theory without particles, which 
is k = 1.53 ±0.10 at /3l = 10 fl25"f . The three results, rows (2), (5), and (6) of the table, differ 
only in the number of particles used, not in the lattice spacing, the particle charge, or the 
manner in which finite (3^ lattice spacing systematics were taken care of. They are grossly in 
conflict if k is independent of m 2 -,, but they agree very nicely with the ASY scaling law. We 
illustrate this in Figure [F[ Combining them, we get an estimate for k' of k' = 53 ± 3 sta t ± 5 syst . 
We should also fold in the systematic errors from finite lattice spacing and finite Q 2 , which 
we would estimate based on the other runs to be in the 20% range. These then dominate our 
uncertainties, and our final answer is k' = 53 ± 11. For comparison, we can take the results 
for pure Yang-Mills theory without particles from [f25| and extrapolate them to zero lattice 
spacing assuming the ASY scaling applies. The extrapolation is k = 23.6//3l, see figure 7 
of that paper. We can convert this into a value for k' by using the discussion after Eq. 
([Tip, that is, using Arnold's calculation of the relation between damping from hard classical 
lattice modes and from correct hard thermal loops E5|. We get n' — 51 ± 15.^ The error 



here is almost all systematic, arising from the rotational noninvariance of the spectrum of 
lattice modes. 

It is encouraging that these results agree within (admittedly substantial and mainly 
8 We thank Peter Arnold and Jan Smit for discussions on this point. 
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Figure 7: Results for k at three values of particle number but a common value of lattice 
spacing, plotted against g A T 2 /mp. The error bars in reflect uncertainty in the damping 
from hard classical lattice modes. The "old picture" predicts a flat line, while the ASY 
picture predicts a straight line through the origin, like the illustrated fit. 

systematic) errors. It appears that a consistent picture for the Acs diffusion constant has 
emerged. 

8 Conclusions 

We have developed a procedure for generating the nonlocal "hard thermal loop" effects in 
classical field simulations of Yang-Mills theory, by introducing particle degrees of freedom 
which generate these effects but can be treated with reasonable numerical effort. The new 
system is Hamiltonian and has the same thermodynamics as thermal Yang-Mills theory in the 
dimensional reduction approximation. We have also tested that the self-energy corrections 
due to the particles in the abelian theory are correct, and we have analytically computed 
the plasma frequency at leading order and it is also correct. 

Further, we have applied the technique to the calculation of the diffusion constant of 
Chern-Simons number in pure Yang-Mills theory. Our results vindicate the proposed scal- 
ing law of Arnold, Son, and Yaffe, that the diffusion constant T should scale inversely with 
the square of the Debye mass. This statement is not sensitive to the quality of the ther- 
modynamic corrections we have applied, because we have data with different numbers of 
particles which show starkly different diffusion rates, but the lattice spacing and hence the 
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thermodynamic corrections we have applied are the same. Our results (in Yang-Mills theory, 
which will be valid in Yang-Mills Higgs theory only for thermal Higgs mass m%(T) ^> g 4 T 2 ) 
convert to T = 29 ± 6« 5 T 4 at ml = llg 2 T 2 /Q. 

There are no obstacles to extending our techniques to Yang-Mills Higgs theory or groups 
larger than SU(2), which would allow us to determine the screening dependence of a number 
of other dynamical properties of interest for baryogenesis, such as the IR bosonic contribution 
to the bubble wall friction and the strong sphaleron rate. The technique may also have 
applications for relativistic heavy ion collisions. 
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A Plasma frequency in the abelian theory 

We have shown that the thermodynamic modifications due to particles are very close to 
those expected, and we have argued that the role of the particles in generating hard thermal 
loops should be the same as in the continuum theory, for suitably infrared and slowly varying 
gauge fields. Here we will check this, and test its limits, by explicitly calculating the plasma 
frequency, the oscillation frequency of a spatially homogeneous electric field, in the abelian 
theory. The calculation should apply approximately to the nonabelian theory, in the regime 
where u p ^> g 2 T, since in this case the theory looks approximately abelian on the relevant 
length scales. Away from this limit, the plasma frequency of the nonabelian theory is not 
well defined, since the electric fields only oscillate coherently on a time scale of order l/g 2 T. 
We will calculate in the formal small Q, large (n), but finite Q 2 (n) limit, but we will deal 
explicitly with the discrete nature of the lattice electric fields and their interactions with 
particles. 

First we will consider the case of a continuous electric field. It is sufficient to consider 
an electric field pointing along a lattice direction, since we will work to linear order in the 
field strength, so general fields can be studied as a linear combination. We take the field 
to be Ei(t) = 5 ix ReE e luJt , with E a constant giving the strength of the electric field, and 
we will solve self-consistently for u. Consider a particle of charge q <C 1 propagating in the 
background of this field, with a mean momentum in the x direction of p x and a momentum 
orthogonal to the x direction of p±. The x momentum will have a time dependent disturbance 
of Sp x (t) oc qE , due to the electric field, which satisfies 

5p x = qE x = ReqE e iuJt -> 5p x = Re^^e iuJt . (75) 
The electric field responds to a current in the x direction due to this particle, of magnitude 
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qv x . Expanding in 8p x <C 1, This is 



Px + SPx Px , fax (, Pl\ Px . 5PxP± f-g* 

+ — J~ ■ ( 76 ) 



pl+p 2 x + 2px6px P P V r) P P 



Now averaging over particles, the leading term will cancel between particles of opposite sign 
for q, and the time derivative of the electric field will be 

_ 4 i E s*£ = Re df*,« E d. (77) 

particles ^ " particles " 

Now E x = Keiu E e luJt . The E cancel on the two sides, and we get an expression for uj 2 . A 
nice way of thinking of this is that what we want to know is Ke(—iE/E) = uj; we just plug 
the determined value of E into this expression to find u. 

Now we evaluate the expression. The average over angles of p\ is 2p 2 /3, and the average 
over particles of 1/p is 1/2T. 

2 = \ /W / 7g 
3T v ; 

In lattice units, 1/T becomes In the nonabelian case, (g 2 ) = Q 2 /3, because that is the 
average of the square of the projection of a particle charge into one abelian subgroup. We 
therefore find uj 2 = 7^/3, as it should. 

Now let us include the discrete nature of the electric fields and their interactions with 
the particles. Again we take the electric field to be spatially uniform with the same value, 
and we will solve self-consistently for the plasma frequency. 

Again, what we want to know is Ke(—iE/E), or its time average. We consider the 
influence of one particle, then we will average over all particles. Again the particle has a 
mean momentum in the x direction of p x and a perpendicular momentum of p±. We will 
often write ± to mean sign^) and =F to mean —sign(p x ). 

The impulse on the electric field on some link, due to a particle crossing the dual lattice 
face the link penetrates, is of fixed modulus =Fg. Past interactions with E change the current 
the particle induces only by changing the time at which the impulse is felt .0 Let the particle 
passage time, at zero order in q, be t CTO ss- This crossing time will receive a correction due 
to past interactions of the particle with the electric field, which we denote St(t CTOSS ). The 
particle's contribution to Re / dt(—iE/E) is 

E Re J dt(-iE / E) = Re[(±iq)e- iuJ ^° ss+St(t " oss)) ] (79) 



Re 



(±iq)e~ luJtcioss +Re ±qu6t(t croBB )e 



-'LUjtr 



The leading term is q odd, and cancels between particles of opposite charge. But St(t cros8 ) 
will be q odd at leading order, and this term will contribute at 0(q 2 ). 

9 Thcre will also be a small number of particles which would have induced a current on this link, but 
instead induce a current on its neighbor. But there are an equal number which induce a current on this link 
rather than its neighbor, and this cancels out. 
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Let us compute 5t(t) to leading order in q. It is related to the difference in the particle's 
position from the E = case, Sx(t), by St(t) = — 5x(t)\p\/p x . This in turn is given by an 
integral over the particle's past history of the difference in its velocity from the average value, 



5x(t)= f dt'5v x (t') = f dt'5p x (t')^. (80) 

J — CO J —oo p 



Now 5p x (t') is a sum over past crossings of the kick to the particle at that crossing. The 
particle has undergone an infinite number of past crossings, one for each plane parallel the 
plane of the dual lattice face it is currently crossing. The time since it crossed the plane a 
distance =Fm away in the x direction is ±m\p\/p x , and the phase of the electric field at that 
time was exp(iu(t CTOSS =p m\p\/p x )), and the impulse it received was, at leading order in q, 
±qE\p\/p x . Hence, 

Across)) = Re (e**™ £ ^big^obi pl e ^ mN/p \ (gl) 

V m=l ^-Vx ±Px P J 

This sum is not absolutely convergent. So to make sense of the calculation, we must 
regulate it. We do this by assuming that, as well as fluctuating with the plasma frequency, 
the correlator of E at two times has a slow exponential envelope, so E x in the past is not 
Re£ e^ but ReE e iiu)+t)t . Of course, we will take the e -> limit at the end. 

Thus, after regulating the sum, the kick from this particle crossing this boundary is 



Re I dtl^L = Re f; T !H^^ e T^) m \ P \/ Px (g2) 

J E m=l Px 



The frequency with which a particle of this momentum encounters a boundary in the x 
direction is ±p x /p, so u is 



,2, ,J2 



- J2 ReJ^-^^e^^^. (83) 

^ particles m=l Px\P\ 



Now writing x = p x /\p\, and performing the integral over p 2 dp and over the aximuthal angle, 
this reduces to 



uj 2 = ~ ^ '^ Relim 

3 e-»o 2 



The sum over m may now be performed, 



/ dx— -=- V me { - tuJ - e)m/x . (84) 
Jo x 2 ^ } 



£ me^ 1 = - sinh- 2 ■ (85) 



m=l 

Making the substitution y — 1/x, the integral becomes 



uj 2 = x "' ^ /r ^ u Re lim 



6-o 2 A y Ay 2 \ 2 j 
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For e > the integral is exponentially convergent, and we can rotate the contour to run in 
the negative imaginary direction: 



2 = \ /vi '^ Relim / idy y \ ' ' sinh^ ^ ^ i . (87) 

3 2 Jo 4(?/ + f) 2 \ 2 J K ' 

The integral is well behaved and we are free to take the e — > limit. Taking the real part, 
after some algebra we obtain 

^2 = MM^ xF( ^), (88) 
3 

F( 2 \ — f°° a %( cosn (^2/) cos(^) - 1) + (?/ 4 + 3y 2 ) swhjujy) sin(o;) 

[UJ ' ~ 2 Jo V 2(l + 2/ 2 ) 2 (cosh(a; ? /) - cos(u;)) 2 1 J 

w 2 u; 4 a; 6 
^ 1 - -r + ^tt: + 



4 240 30240 

Here, of course, u is the plasma frequency in lattice units, ie auj in physical units. We see 
that it receives a correction when u ~ 1/a, but that for u; 2 a 2 ~ 1/2 in lattice units, the 
correction is already quite small. There will also be 0((n)Q 4 ) corrections, probably including 
damping, but we have not been able to calculate these analytically. 



B Integrating out the Aq field at arbitrary mo 

We would like to understand the relation between the thermodynamics of infrared magnetic 
fields in our lattice system and the thermodynamics of infrared magnetic fields in the full 
quantum theory, say in the dimensional reduction approximation. As we have seen, our 
system has thermodynamics described by a lattice gauge theory with an A field, which has 
a Debye mass which depends on /3l and on {n)Q 2 , whereas the dimensional reduction 
approximation gives the continuum limit of this lattice theory, with some particular renor- 
malized determined by the particle content of the full quantum field theory. There are 
two complications here, the difference in behavior between lattice and continuum systems 



and the difference in Debye mass. The first has been dealt with in pT| , for the system 
without the A field and for the system with the A field but in the approximation that mp L 
is small in lattice units. In the theory with only classical lattice modes this approximation 
is parametrically justified, but with the inclusion of particles this is not necessarily the case. 

To deal with finite m^, one first notes that is large enough to make the Aq field 
heavy, and since the theory is super-renormalizable one can integrate over such heavy fields 
at one loop and capture their dominant contributions to the infrared. We should perform 
this integration in each system, leaving us to compare pure Yang-Mills theory in 3-D, on the 
lattice and in the continuum. The matching between these has been studied in [[HJ up to 
corrections of order fi^ 2 , so here we only discuss the integration over the Aq field. 

Integrating over the A field in the continuum theory was studied in || Bl . There is only 
one correction, a self-energy correction to the gauge fields which shifts the gauge coupling 
gl = g 2 T to 
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(a) (b) (c) (d) 

Figure 8: Diagrams needed in the integration over the A field in Yang-Mills Higgs theory. 
The double solid line is the A field, the single solid line is the Higgs field, and the wavy line 
is the gauge field. 



In the theory with a Higgs field there is also a shift to the scalar self-coupling, 

A 3 = A 3 - T7 r|^, (91) 
1287rmD 

and to the Higgs mass squared parameter, 

-3 = ™ 2 - ^ + 0( 9 t/16^) . (92) 

The scalar wave function is not corrected. 

In the lattice theory the integrals are trickier because the A field has lattice dispersion 
relations and the gauge- t4 vertices have nontrivial momentum dependence. So to warm up 
we will start with the corrections to the Higgs parameters. Our unit conventions will be the 
same as in [11], ie the lattice scalar self-coupling will be Al = 4X/g 2 , and the lattice spacing 
will appear in the dimensionless quantity /3l = 4/(g 2 aT). 

Two diagrams are relevant, Fig. || diagrams (a) and (b). There are also high dimension 
operators induced by diagrams with more <fi lines, but the influence of these terms is small 
for large m 2 ^ and we can drop them, just as in the continuum case ||. Diagram (a) gives 
a momentum independent self-energy correction, ie a mass squared correction. Denote the 
coupling between the A field and the Higgs field as Aa; at lowest order Aal = 2. Here and 
in what follows an L subscript means the value in lattice units, with the normalizations used 
TTJ. The shift in the lattice Higgs mass is^| 



in 



_ 2 2 3/iAL f d 3 k 



2(j L J[-7r,7r]^ (Any k 2 + mf, L 
where hi = 2sin(A;j/2) and 

k 2 = J2tf = XX 2 -2cosJfei). (94) 

i i 

Henceforth the range of integration and the d 3 k/(2^) 3 will be understood. 

10 Note that there are several typos in the Appendix C of in which factors of S/47T or £/47r are left 
out. 
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The integral in (|93|) must be determined numerically, and we define it as 



1 _ E(m DL ) £ m DL £m DL m DL 

+ + ■ ■ • , ( 95 ) 



fc A; 2 + m DL 4tt 47r 47r 47r 327T 

where £ = S(0) = 3.175911536 and £ = 0.152859325. We display the small ttt-dl expansion, 
but at general m^L the integral should be done numerically. 

At zero external momentum diagram (6) corrects the scalar self-coupling by 

* l=Al -i& /* (^+Ul) 2 ' (96) 

We define 



1 = ^dl) ^ 1 + J_ _ 3m DL + 

/fc (/c 2 + m DL ) 2 47r 87rmDL 47r 647r 

Again, we have displayed the leading terms at small mDL, but at general moL the integral 
should be treated numerically. This definition of £ (ttjdl) does not have £(m — > 0) — > £; 
instead £(m) has a pole at m = 0, and £ is the constant part after the pole. The residue of 
the pole reproduces the continuum expression, (pTj) , and the presence of the pole is because 
we are computing the full effect of the A field on the lattice, not the difference between 
lattice and continuum theories. 

Now we treat gauge field self-energy corrections. Since there are no A field corrections 
to the three point gauge-scalar vertex or to the Higgs wave function at one loop, this is the 
only further correction we need. There are two diagrams, (c) and (d) of Figure || Denoting 
the polarization and spin indicies as i,j and a, b, at general external momentum p, diagram 
c gives 



8S ab Sij f cos hi 



Pl J k 2 + m DL ' 



and diagram (d) gives 



S a b f 16 sin ki sin kj 



PlJ [(k-p/2)* + ml h }[(k + p/2y + ml h } 



{9i 



(99) 



The sum of the contributions at p = is 

85ij cos ki 16 sin ki sin kj \ 5 a b f 9 




-'tlim^^^- (100) 



k 2 + m 2 DL (k 2 + m DL ) 2 J Pl 

The last equality is because the integrand is a total derivative without singularities and the 
domain of integration is compact without boundary. Hence the A field does not induce a 
mass for the gauge field, as ensured by gauge invariance. 

Next, we find the 0(p 2 ) term in the self-energy, which is responsible for renormalizing 
the coupling. Only diagram (d) depends on p, and expanding Q99D to second order gives 



46 



AB 



Ik 



^ / r 2 sin fa sin kt cos ki5i m f 4 sin fa sin fe,- sin fa sin k m \ , 
im \ J (k 2 + m^ L ) 3 J (P + m DL ) 4 J 
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To evaluate this, we need the identity 



Q2 



cos ki 



dhdkj k 2 + 



DL 



8 sin ki sin kj cos k\ 



(k 2 -mj 



25ij cos ki cos ki cos kiSiji 



{k 2 + m 2 -^ 2 



b VL) ~r ""T)L/ k 2 + m^L 

where we have already used ( |100| ) to simplify the final expression, and the identity 
«9 4 



(102) 







dkidkjdktdk^^ +mcJ 
96 sin ki sin kj sin ki sin k m 

(jfc2 + m 2.U 



+ 



]_g^' cos ^* s * n ^ s * n ^ m + ^ permutations 



+ 



+ 



65 ijM cos fcj 



DL; " ' '"'DL 

4 cos^cosA;,^^ + 2 permutations 



DL 



(A; 2 + m 2 ^) 



(103) 



Further, we use the relationship 

cos fcj - 

from which it follows that 

cos h 



DL 



m 



DL 



(104) 



S(m DL ) 



k 2 + 

and, with a little more work, that 
cos ki cos 




(105) 



(//-' • m 2 



b DL) 



1 + IUhl + m DL^ £(™Dl) 



24 



47T 



'I + !^DL^ S ( m DL) 

4 24 



' ^DL + m Dh\ £(™Dl) 



21 



47T 




(106) 



Using these, after considerable algebra we find that the 0(p 2 ) contribution to the self-energy 
is 



Sab (, 2 
PL 



4 _ 2ml L _ m^A £(m DL ) .A, ™dl\ S(m DL ) 
3 3 18 / 4tt 1 3 18 / 4vr 



(107) 



While intermediate expressions have been rotationally non-invariant, the result is rotation- 
ally invariant and transverse. This happens because the only gauge invariant dimension 4 
operator which is cubic invariant is F 2 -. At 0(p A ) there will be rotationally non- invariant 
contributions, although they will be cubic invariant. We are not concerned with these here 
because we are not trying to compute induced nonrenormalizable operators, only corrections 
to terms already in the Hamiltonian. But the nonrenormalizable terms could be computed 
by a straightforward but extremely tedious extension of what we have done here. 
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yL, naive 




2 

Tfl ^ 

'"D 




P { Tnr\ 1 


ML, imp 


8 7 


375 


1 59 


1 901 


3487 

. kj t: O i 


8 073 


10.7 


.189 


1.09 


2.101 


.4556 


10.078 


12.7 


.111 


0.74 


2.280 


.5821 


12.085 


10.7 


.378 


1.77 


1.840 


.3214 


10.072 


10.7 


.756 


3.14 


1.498 


.1979 


10.069 



Table 2: Conversion between /^L.naive and Avmp for the simulations in this paper. Everything 
is in lattice units. In Section [7|we truncate Ph,imp ^° the nearest integer when we write it. 



The rriDL = limit of ( |107| ) agrees with the result in [|TT[, and the coefficient of the pole 
at 777.DL = reproduces (|90|) . In the opposite limit of very large moL one can (Taylor) expand 
£ and £ in m^L, 



E(m) 

47T 

Air 



m 



6m -4 + 42m" b - 324m" 8 + 2610m" 



-10 



m~ 4 - 12m~ b + 126m" 



1296m" 10 + 13050m" iZ - 



12 



(108) 
(109) 



The expansion converges for m 2 > 12 = sup(k 2 ). While the £ and E contributions in (|107| ) 
separately give 0(m°) contributions, these cancel, as do the first several powers of m 2 , and 
the correction is S5 a b(p 2 Sij — ViVj) I {PiJ^h) a ^ leading order in m^L- So the A field rapidly 
becomes irrelevant at large Debye mass, as it should. 

Combining the result from integrating out the A Q field with the thermodynamic correc- 



tions from transverse gauge boson loops found in ||11|| , we get a relation between the "naive r 
tree level /3L,naive and the improved value, 



II 37£\ 
ML, naive PL,imp I „ „ I ^ 

\ 3 67r / 



4 
3 



2ml h 



£(™dl) 



m£> L \ E(m DL ) 




(110) 



Here the constant £ = 0.152859325 is the limit of £(m) as m goes to zero, after the pole has 
been removed. 

For completeness, we also list the values of Ph,imp an d A^naive for the simulations presented 
in section [7], in Table 

This completes the integration over the Aq field at one loop, at general m 2 ^. 



C Probing the retarded propagator by linear response 

In this appendix we show how to probe the retarded photon propagator by studying the 
linear response of the abelain plasma to an external current. 

First we consider the response to transverse perturbations. In the framework of linear 
response theory, the response A(x) is related to the external current j(x) through the retarded 
propagator D R (x,x'): 

Ai(x) = J d A x'D R (x,x') Jk (x') . (Ill) 
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Going to Fourier space and suppressing the Lorentz indices, we have 

E(u, k) = iuD R (uj, k)j(u, k) , 

where 



D R (u,k) = 



1 



k 2 - U(cu, k) ' 



or 



(112) 
(113) 



Taking ji(x,t) = jijoe tlJJot e 9(t) as the driving current, we find 



Ei(u,l) =ji 



3ouj(27r) 3 5 3 (k-l) 



(to — uj + ie)(uj 2 — k 2 — H(u, k)) 
Fourier transforming E(u>,k) back to space-time gives 



(114) 



E(x,t) = j e ikox J 



ooe 



-iurt 



2tt (lu - lu + ie)(uj 2 - k 2 - U(u, k)) 



(115) 



The retarded propagator D R (u, k) is analytic in the upper half complex w-plane due to 
causality, and only has poles in the lower half plane. Writing H(u, k) = U t (uj, k) + iU[(uj, k) 
and noting that H r (—ou, k) = Tl r (ou, k) and IIi(— ou, k) = — Hi(ou, k), the two poles are located 
at u} p i(ko) = uj p (ko) — i'jp(ko) and — u^ko), where u p (ko) is the plasma frequency for the 
mode ko and 7 P (fco) is the associated on-shell damping rate. 

Completing the contour integral in the lower half cu-plane, we find 



E(x,t) = 



joUJ e- iuJot e ik ° x 
w o - ko - n(w , h) 



0(t) 



uJ p i(k )-Uo u^i^+ujo 



6{t) . (116) 



The first term represents the asymptotic response, while the second term is the plasma 
oscillation at k — ko, which is only transient due to damping at a rate of 7 P (fco). At large t 
(t 7p X (A;o)), after the transients have died out, one sees only the asymptotic behavior 



E(x, t) 



. jou} e- iuJot e ik ° x 
' Uo -kl~ n(w , ko) 



117) 



A current of the form j(x,t) = j sin (u; t) sin (k x) • Q(t) will generate the following E 
field at large t, 

E(x, t) = joUJoP cos (uoot + a) sin (k x) ■ 9(t) , (H8) 



where 



p = {[ul-kl-U r (u ,ko)} 2 + U 2 (u ,ko)}' 
-U.i(u ,ko) 



a 



tan 



u 2 - H - Il r (u;o, k ) 



(119) 
(120) 
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The real and imaginary part of H(uj, k) can therefore be determined from the space-time 
averages 

a; ni(a;o, k ) 



and 



(j(x,t)-E(x,t)) 



(j(x,t)-E(x,t-r)) 



jI 



4 K 2 - kl - n r (w , k )} 2 + n 2 (o; , k ) 



1 - kl- II r (a;o, fc )] 



121 



(122) 



4 [ul - kl - U T (u , k )} 2 + U 2 (u , fco) ' 
where r = ir/(2uo), and the averages are taken after the decay of the transients. 

The numerically determined values for the transverse self-energies Ili and Il r can then be 



compared to their perturbative values (48 

3 



U t (uj, k) 



ILi(u,k) 



2^F 



1 + 



l,k 



uj. 



2 k uj 



In 



uj + k 



3tt 2 uj 2 k 



UJ 



W(k 2 



uj — k 
uj 2 ). 



(123) 



(124) 



K u k' 

The response to longitudinal perturbations is easiest to treat in Coulomb gauge, where 
only Ao is nonzero and E is determined from diA$. The retarded propagator is 

,2 ^ 



One obtains 



Df (uj,k) = 

(j(x,t) ■ E(x,t)) 
(j(x,t)-E(x,t-r)) 



or 



k 2 (uj 2 -k 2 -U L (k,u;)) ' 

j 2 {uj 2 - fc 2 )n Lji 



2 ' 



4uj (uj 2 - k 2 - n Ljr ) 2 + n 2 . 
io (uj 2 -k 2 )(uj 2 -k 2 -n Lir 



Auj (uj 2 - k 2 - n L , r ) 2 + n 2 ;i 

The real and imaginary parts of the longitudinal self-energy are 



n L)r (a?,/c) 
n L ,i(o;, k) 



3^(1 



1 



UJ 

2k 



In 



uj — k 



3n 2 uj uj , 2 2 



(125) 

(126) 
(127) 

(128) 
(129) 



D Relating time scales when there are particles 



It has been shown |jTT|, as discussed in Appendix [FJ, that the difference in screening from 
UV modes between lattice gauge theory and continuum gauge theory can be understood at 
leading order in a (or in P^ 1 ) as a rescaling of the differential operator D i: so what at tree 
level looks like D 2 A in fact behaves like Z^D 2 A, with Za = 1 + 0(a). The 0(a) term has 
been computed in |ll] and extended to arbitrary Debye mass in Appendix [FJ. Hence, in the 
case with no particles, the equations of motion for Ai in temporal gauge look like 

dE- 1 

' = -Z2 1 (D%--(D i D j + D j D i ))A j , (130) 



dt 
OA 
dt 



ZeE% 



(131) 
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and (E 2 ) oc /?l n aive . The usual rule A = E is rescaled by which has not been computed, 
but which lacks the large "tadpole" contributions which characterize Za- To get a complete 
0(a) correction of the dynamics it would be necessary to compute Ze, but it has been 
advocated in [^] that the absence of "tadpole" contributions means that Ze — 1 can be 



1/2 

neglected compared to Za — 1. In this case, a simple rescaling of time, Z A d/dt = d/dt', 

1/2 

and of the electric field, Z/ E = E , eliminates Za, and replaces /3L,naivc with /3^i mp = 
■^A 1 /^L,naive- So the time scale used, t, is related to the more appropriate time scale, t', by 
t' = tw/5L,imp//3L,naive- When one computes T one should divide (Nq S ) by Vt' rather than 



Vt, and the corrected rate is J /^L.naive/ Plmp times larger. 

However we are now interested in the case where there are enough particles to put the 
evolution in the overdamped regime. Huet and Son [p4] and Son ]|9] have recently shown 



what new term this adds to ( |130| ). A thermal distribution of particles would not contribute, 
since particles of opposite charge would cancel, but the particle distribution at any point is 
skewed because of past, remote electric fields, and contains fluctuations. The new evolution 
equation looks like 

+ Z A \D% - \{D iDj + DjDi))Aj(x,t) = 9M^1 

x J dyN$(x, y, t)E)(y, t - \x - y\) + (noise term) . (132) 

Here the right hand side represents the departure of the particle population from rotational 
invariance, due to its linear response to past electric fields, and fluctuations in the particle 
population. The form of the nonlocal kernel iV is given in |^4[ fi9|| . They point out that the 



overdamped limit corresponds to the large Q 2 {n) limit, in which case the dE/dt term and 
the time dependence of E on the right hand side can be neglected. Again approximating 
that A = E, rescaling time to eliminate Za and to replace all appearances of /^L.naive with 
Avmp now requires t' = Z~^t. Hence one should compute everything using /?L,im P and the 
lattice spacetime volume used, but then correct the rate by a factor of (/^naive/Avmp) 2 - 
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